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Abstract. We approximate the price of the American put for jump diffusions by a sequence of functions, which 
are computed iteratively. This sequence converges to the price function uniformly and exponentially fast. Each 
qq ■ element of the approximating sequence solves an optimal stopping problem for geometric Brownian motion, and can 

f~ be numerically computed using the classical finite difference methods. We prove the convergence of this numerical 
C ~ ) i scheme and present examples to illustrate its performance. 
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Q , Based on the results of [6] we propose a numerical algorithm to calculate the American option prices for jump 
diffusion models and analyze the convergence behavior of this algorithm. 
' As observed by [6], we can construct an increasing sequence of functions, which are value functions of optimal 

stopping problems (see (2.8) and also (2.11)), that converge to the price function of the American put option 
. uniformly and exponentially fast. Because each element of this sequence solves an optimal stopping problem it 
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1. Introduction 

Jump diffusion models arc heavily used in modeling stock prices since they can capture the excess kurtosis and 
skewness of the stock price returns, and they can produce the smile in the implied volatility curve (see [10]). Two 
well-known examples of these models are i) the model of [18], in which the jump sizes are log-normally distributed, 
and ii) the model of [17], in which the logarithm of jump sizes have the so called double exponential distribution. 



shares the same regularity properties, such as convexity and smoothness, with the original price function. Even the 
corresponding free boundaries have the same smoothness properties (when they have a discontinuity, which can 



■ only happen at maturity, the magnitude of the discontinuity is the same; see (2.13)). Therefore, the elements in 
this approximating sequence provide a good imitation to the value function besides being close to it numerically 
(see Remark 2.1). On the other hand, each of these functions can be represented as classical solutions of free 
boundary problems (see (2.9)) for geometric Brownian motion, and therefore can be implemented using classical 
finite difference methods. We build an iterative numerical algorithm based on discretizing these free boundary 

■ problems (see (3.10)). When the mesh sizes are fixed, we show that the iterative sequence we constructed is 
monotonous and converges uniformly and exponentially fast (see Proposition 3.3). We also show, in a rather direct 
way, that when the mesh sizes go to zero our algorithm converges to the true price function (see Proposition 3.4). 

The pricing in the context of jump models is difficult since the prices of options satisfy integro-partial differential 
equations (integro-pdes), i.e. they have non-local integral terms, and the usual finite-difference methods are not 
directly applicable because the integral term leads to full matrices. Recently there has been a lot of interest in 
developing numerical algorithms for pricing in jump models, see e.g. [1], [2], [4], [9], [11], [13], [14], [16], [17], [19], 
[23], among them [2], [9], [13] and [14] treated specific or general jump models with infinite activity jumps. These 
algorithms have been extensively discussed in Chapter 12 of [10]. In this paper, relying on the results of [6] as 
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described above, we give an efficient numerical algorithm (and analyze its error versus accuracy characteristics) to 
efficiently compute American option prices for jump diffusion models with hnite activity. One can handle infinite 
activity models by increasing the volatility coefficient appropriately as suggested on p. 417 of [10]. 

An ideal numerical algorithm, which is most often an iterative scheme, *should monotonically converge to the 
true price uniformly (across time and space) and exponentially fast*, that is, the error bounds should be very 
tight. This is the only way one can be sure that the price output of the algorithm is close to the true price 
after a reasonable amount of runtime and without having to compare the price obtained from the algorithm to 
other algorithms' output. It is also desirable to obtain a scheme **that does not deviate from the numerical 
pricing schemes, such as finite difference methods, that were developed for models that do not account for jumps**. 
Financial engineers working in the industry are already familiar with finite difference schemes such as projected 
successive over relaxation, PSOR, (see e.g. [20]) and Brennan-Schwartz algorithm (see [8] and [15]) to solve the 
partial differential equations associated with free boundary problems, but may not be familiar with the intricacies 
involved in solving integro-partial differential equations developed in the literature. It would be ideal for them if 
they could use what they already know with only a slight modification to solve for the prices in a jump diffusion 
model. In this paper, we develop an algorithm which establishes both * and **. In Section 4, we will name 
this algorithm, depending on which classical method we use to solve the sparse linear systems in (3.10), as either 
"Iterated PSOR" or "Iterated Brennan-Schwartz" . 

In the next section we introduce a sequence of optimal stopping problems that approximate the price function of 
the American options, and discuss their properties. In Section 3, we introduce a numerical algorithm and analyze 
its convergence properties. In the last section we give numerical examples to illustrate the competitiveness of our 
algorithm and price American, Barrier and European options for the models of [17] and [18]. 

2. A Sequence of Optimal Stopping Problems for Geometric Brownian Motion Approximating 



We will consider a jump diffusion model for the stock price St with Sq = S, and assume that return process 
X t := log(St/ S), under the risk neutral measure, is given by 



In (2.1), fi = r + A — A£, r is the risk- free rate, Wt is a Brownian motion, Nt is a Poisson process with rate A 
independent of the Brownian motion, Zi are independent and identically distributed, and come from a common 
distribution F on R, that satisfies £ := J R e z F(dz) < oo. The last condition guarantees that the stock prices have 
finite expectation. We will assume that the volatility a is strictly positive. The price function of the American put 
with strike price K is 



in which St,T is the set of stopping times of the filtration generated by X that belong to the interval [t, T] (t it the 
current time, T is the maturity of the option). Instead of working with the pricing function V directly, which is 
the unique classical solution of the following integro-differential free boundary problem (see Theorem 3.1 of [6]) 



the American Option Price For Jump Diffusions 




(2.1) 



V{S,t):= sup E{e- r(T -*)(A"- S T ) + \S t = S}, 



(2.2) 




V{S,t)=K-S, S<s(t), 
V(S,T) = (K~S) + , 



(2.3) 
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in which, A is the differential operator 

a -=?' s *w* + » s t s - < 2 - 4 > 

and t — > s(t), t G [0, T], is the exercise boundary that needs to be determined along with the pricing function V; 
we will construct a sequence of pricing problems for the geometric Brownian motion 

dS° = nS°dt + aS°dW t , S§ = S. (2.5) 

To this end, let us introduce a functional operator J, whose action on a test function / : R + x [0, T] — > K + is the 
solution of the following pricing problem for the geometric Brownian motion: (S®)t>o 

Jf(S,t)= supe//" e-( r+A ^ u - t U-P/(SS > «)d« + e~ (r+A)(T ~ t) (Jir-^) + |^ = s}, (2.6) 

reS t ,T J 

in which 

Pf(S,u) = [ f(e z ■S,u)F(dz)=E[f(e z S,u)}, S>0, (2.7) 



for a random variable Z whose distribution is F, and St : r is the set of stopping times of the filtration generated 
by W that take values in [t, T]. Let us define a sequence of pricing functions by 

v (S,t) = (K - S) + , v n+ i(S,t) = Jv n (S,t), n > 0, for all (S,t) G R+ x [0,T]. (2.8) 

For each n > 1, the pricing function i; n is the unique solution of the classical free-boundary problem (instead of a 
free boundary problem with an integro-diffential equation) 

JU„(S,i) + Av n {S,t) - (r + A) • «„(5,t) = -A • (P« n _i)(5,t), 5 > s n (t), 
v n (S,t)=K-S, S<s n (t), ( 2 - 9 ) 



n (S t T) = (K-S)+, 



in which t — > s n (f) is the free-boundary (the optimal exercise boundary) which needs to be determined (see 
Lemma 3.5 of [6]). Now starting from vq, we can calculate {v n } n >o sequentially. For v n , the solution of (2.9) 
can be determined using a classical finite difference method (we use the Crank-Nicolson discretization along with 
Bernnan-Schwartz algorithm or PSOR in the the following sections) given that the function v n -\ is available. The 
term on the right-hand-side of (2.9) can be computed either using Monte-Carlo or a numerical integrator (we use 
the numerical integration with the Fast Fourier Transformation (FFT) in our examples). Iterating the solution for 
(2.9) a few times we are able to obtain the American option price V accurately since the sequence of functions 
{"n}n>o converges to V uniformly and exponentially fast: 

v n (S,t)<V(S,t)<v n (S,t)+K(l-e-^ + ^ T -^y (-^ , SGR+, te (0,T), (2.10) 

see Remark 3.3 of [6]. Note that the usual values of T for the traded options is 0.25, 0.5, 0.75, 1 year. 

Remark 2.1. The approximating sequence {v n } n >o goes beyond approximating the value function V. Each v n and 
its corresponding free boundary have the same regularity properties which V and its corresponding free boundary 
have. In a sense, for large enough n, v n provides a good imitation ofV. Below we list these properties: 
1) The function v n can be written as the value function of an optimal stopping problem: 

v n (S,t):= sup E{e-^ TA ^- t \K-S TA<Tn )+\S t = S}, (2.11) 
in which <j n is the n-th jump time of the Poisson process N t . 
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2) Each v n is a convex function in the S-variable, which is a property that is also shared by V . Moreover, 
the sequence {v n } n >o is a monotone increasing sequence converging to the value function V (see (2.10)/ 

3) The free boundaries s(t) and s n (t) have the same regularity properties (see [7] and its references) : 

a) They are strictly decreasing. 

b) They may exhibit discontinuity at T : If the parameters satisfy 

r < A f (e z - l)F(dz), (2.12) 
Jr + 



we have 



lim s(t) = lim sjt) = S* < K, n > 1, (2.13) 

t->T t->T 



where S* is the unique solution of the following integral equation 



rK + A 



(K - Se z ) + -{K- Se z ) F{dz) = 0. (2.14) 



We will see such an example in Section 4, where the equation (2.14) can be solved analytically for 
some jump distribution F. 
c) Both s(t) and s n (t) are continuously differ entiable on [0, T). 

3. A Numerical Algorithm and its convergence analysis 

3.1. The numerical algorithm. In this section, we will discretize the algorithm introduced in the last section 
and give more details. For the convenience of the numerical calculation, we will first change the variable: x = log 5*, 
x(t) = logs(£) and u(x,t) = V(S,t). u satisfies the following intcgro-differential free boundary problem 

d Id 2 ( 1 n\ d 

—U + 2 CT2 ^2" + ( M- ^ j -g^ U - ( r + X) U + ^ ■ i Iu )( X ,t) = 0, X > X(t) 

u{x,t)=K-e x , x<x(t) ( 3 - 1 ) 
u( Xl T) = (K-e*) + , 



(Iu)(x,t) = u{x + z,t)p(z)dz, (3.2) 



in which 



with p(z) as the density of the distribution F. Similarly, u n (x,t) = v n (S,t) satisfies the similar free boundary 
problem where u in (3.1) is replaced by u n in differential parts and by u n -i in the integral part. In addition, it 
was shown in Theorem 4.2 of [21] that the free boundary problem (3.1) is equivalent to the following variational 
inequality 

C D u{x,t) + A • (Iu)(x,t) < 

u(x,t)>g(x) (3.3) 
[C D u{x, t) + A • {Iu)(x, t)] ■ [u{x, t) - g(x)] = 0, (x, i)elx [0, T], 

in which 

, 9 i 2 d 2 ( i 2 \ a 
Cdu =m u+ 2 a dx^ u+ ^-2 a )d^ u - {r + x)u 

g(x) = (K - e x ) + . 
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Since the second spacial derivative of u does not exist along the free boundary x{t), the variational inequality (3.3) 
does not have a classical solution. However, Theorem 3.2 of [21] showed that u is the solution of (3.3) in the 
Sobolev sense. In the same sense, u n {x,t) satisfies a similar variational inequality 

C D u n (x, t) + A • (Iu n -i)(x, t) < 

u n (x,t)>g(x) (3.4) 

[C D u n {x, t) + A ■ {Iu n _{){x, t)\ ■ [u n (x, t) - g(x)} = 0, {x, ijelx [0, T]. 

Let us discrctize (3.3) using Crank-Nicolson scheme. For fixed At, Ax, x m i n and x max , let M At = T and 
LAx = x max — x m i n . Let us denote xi = x m i n + lAx, I = 0, • ■ ■ , L. By u ' m wc will denote the solution of the 
following difference equation 

- 0p_vf- 1,m + (i + e Po )u l ' m - e P+ u l+1 ' m - b l ' m > o 

u l ' m > 9 l (3.5) 



-ep-u 1 - 1 '" 1 + (1 + 6p )u l > m - P+ u l+1 - m - b l ' m 



[W -g 1=0, 



for m = M — 1, • • • , 0, / = 0, • • • , L, satisfying the terminal condition u LAI = g l = (K — e x ') + and Dirichlet 
boundary conditions. 8 is the weight factor. When 8 = 1, the scheme (3.5) is the completely implicit Euler scheme; 
when 9 = 1/2, it is the classical Crank-Nicolson scheme. The coefficients p_, p + and po are given by 



1 2 At If I 2 \ At 



p - = r jAxy-^^-r j ax> 

1 2 At If 1 2 \ At (3.6) 
= 2 JKxY 2 V 2 J Ax' 
p Q =P-+p+ + (r + A)At. 

The term b is defined by 

gi, m = ^ _ 6)) p _ui-bT»+i + 9)^)2'."*+! + (l - 6»)p + w' +1 ' m+1 + AAt ■ [(1 - 6>)(Iu)'' m+1 + 8(Iu) l > m ] . (3.7) 

/ in (3.7) is the discrete version of the convolution operator / in (3.2). It will be convenient to approximate this 
convolution integral using Fast Fourier Transformation (FFT). Discrctizing a sufficiently large interval [z m i n , z ma x] 
into J sub-intervals. For the convenience of the FFT, we will choose these J sub-intervals equally spaced, such 
that JAz = z max — z m i n . We also choose Ax = aAz, where a is a positive integer, so that the numerical integral 
may have finer grid than the grid in x. Let zj = z m i n + jAz, j = 0, • • • , J. I is defined by 

(iuj = Uinterp (%l + z j , TO At) p(Zj)Az, (3.8) 

3=0 

in which the value of Uinterp is determined by the linear interpolation u. That is if there is some V satisfying 

xv < xi + zj < xi' + i, 

then 

Uinterp (x t + Zj,mAt) = (1 - w)u V '™ + WU l ' +1 ' m , 

for some w € [0,1]. On the other hand, if xi + Zj is outside the interval [x m in, x m ax], the value of Uinterp is 
determined by the boundary conditions. Moreover, in (3.8) we also assume 

,7-1 

p(zj)>0, for all j, and ^p{Zj)<l. (3.9) 

3=0 
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Now (3.8) can be calculated using FFT. See Section 6.1 in [2] for implementation details. 

Note that numerically solving the system (3.5) is difficult due to the contribution of the integral term Iu. 
Therefore, following the results in Section 2, we will discrctizc (3.4) recursively (using the Crank-Nicoslon scheme) 
to obtain the sequence {u n }n>o recursively. Let Hq" 1 = g l . For n > 1, u n is defined recursively by 



(1 + 9p )u l n m - 9p + u n 



(1 + 9p )u 



~M-l,m 



> 



(3.10) 



with the terminal condition {t„ = g and Dirichlet boundary conditions. Similar to (3.7), b n is defined by 



i i r =(i - B) P -u 1 - 1 ^ 1 + (i - (i - e) P0 )u\r +1 + (i - 8) P+ u i n +i > m+i 

+ XAt ■ \(l - 0){lu n _ 1 ) l ' m + 1 + 9{Iu n ^) 1 ^ 



(3.H) 



For each n, we will solve the sparse linear system of equations (3.10) using the projected PSOR method (see eg. 
[20]). 

Remark 3.1. We will iterate (3.10) to approximate the solution of (3.5), which can be seen as a global fixed point 
iteration algorithm. This global fixed point algorithm is different from the local fixed point algorithm in [11], where 
d'Halluin et al. implemented the Crank Nicolson time stepping of a non-linear integro-partial differential equation 
coming from an alternative representation (due to the penalty method) of the American option price function. Also 
see [12] for the case of European options. Note that discretizing the non-linear PDE that arises from the penalized 
formulation introduces an extra error. We work with the variational formulation directly. 

Each u n approximates u n , which itself is the value function of an optimal stopping problem, and as we have 
discussed in Remark 2.1 provides a good imitation of the American option price function. Each of these iterations 
provide strictly decreasing free boundary curves with the same regularity and jump properties as the free boundary 
curve for the American option price function, see Remarks 2.1 and 4-. 2. The approximating sequence in [11] does 
not carry the same meaning, it is a technical step to carry out the Crank Nicolson time stepping of their non-linear 
integro-PDE. 

3.2. Convergence of the Numerical Algorithm. In the following, we will show the convergence of the numerical 
algorithm for the completely implicit Euler scheme (8 = 1). We first show that {M n }n>o is a monotone increasing 
sequence. Extra care has to be given to make the approximating sequence monotone in the penalty formulation of 
[11] (see Remark 4.3 on page 341), but the monotonicity comes out naturally in our formulation. Next, we prove 
that the sequence {u n }n>o is uniformly bounded above by the strike price K and converges to u at an exponential 
rate. At last, we will argue that as the mesh sizes Aa; and At go to zero u converges to the American option value 
function u. In the following four propositions, we let At and Aa; to be sufficiently small so that constants P — and 
p + defined in (3.6) are positive. 



Proposition 3.1. The sequence {u n } n >o is monotone increasing sequence. 
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Proof. When 6 = 1, subtracting the third equality for n-th iteration in (3.10) from the equality for (n + l)-th 
iteration, we obtain 



■p+u^-fy 



~L,m 
U n+1 



p_fi}r 1,TO + (i+po)fii 
+ - filr 1,ro ) + a +po) - -*>+ 



~l+l,m _ f,l+l,m 



(3.12) 



u n+l 



= 0. 



in which we used the linearity of the operator I. Let us define the vectors 



e n+l 

era 
Jn+1 



~0,m _ ~0,m . . . -L.m _ ~L,m 



(<#™ +1 - ^™ +1 ) + AAi • ( I(u n - u n _r) 



0,m 



(^ +l -^- m+i )+ 



Equation (3.12) can be represented as 
in which the matrix ^4's entries are 



AAi • (l(u n - u n . 

A m _ fm 



~0,rn 
U n+1 - .9 



-L.m L 
U n+1 - 9 



(3.13) 



at 



(~l,m 
K+i 



(i+po) - + (-p-vir 1 ^ + (i+PaWr -P+u l +^ m - i\r) 



, ~l,m 1 
-P+ [ U n+1 -9 





3 = 1-1 
3=1 
j =1 + 1 
others. 



On the other hand, using the first and second inequalities in (3.10) and the fact that p— and p+ are positive, we 
see that A is an M-matrix, i.e. A has positive diagonals, non-positive off-diagonals and the row sums are positive. 
As a result all entries of A^ 1 arc nonncgative. 

Now we can prove the proposition by induction. Note that u\ > u = g, as a result of the second inequality in 
(3.10) and the definition of uq. Assuming u n > tt n _i, we will show that u n +\ > u n , i.e. u n +i — u l n m > for all I 
and m, in the following. 

. , . / . \ l.m 

is nonnegative from 



First, the terminal condition of u n gives us — u l n M = 0. Second 



\I{u n - u„_i) J 



the assumption (3.9). Assuming ujj™^ 1 — u^ m+1 nonnegative, we have in (3.13) as a nonnegative vector. 

Combining with the fact that all entries of A -1 arc nonnegative, the nonnegativity of u,^ 1 — u 1 ^ 1 follows from 
multiplying A^ 1 on both sides of (3.13). Then the result follows from an induction m. □ 

Proposition 3.2. {ii n }„>o are uniformly bounded above by the strike price K. 

Proof. When 6 = 1, in the third equality of (3.10), there are some (l,m) such that = g. Otherwise we have 

(1 + Po )u l ,r = p-u l ~ l > m + P +u l + hm + «n m+1 + AAi (ln„_i 
However, in both cases, we obtain the following inequality 

(1 + p ) \vk m \ < p-B™ + p+B™ + Bl L+1 + AAtB„_i + rAtK, < I < L,0 < m < M - 1, (3.14) 
in which we define 

B% = (max | ull m | ^ \f A, B n = max B™ . 
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Note that the right hand side of (3.14) is independent of I. Moreover, (1 + po)K is also less than or equal to the 
right hand side of (3.14). Therefore, (3.14) gives us 

(1 + (r + X)At) B™ < B™ +1 + XAtB^ + rAtK. (3.15) 

Given B™ +1 < K and B n _x < K, it clear from (3.15) that -B™ 1 < K. Now the proposition follows from double 
induction on to and n with initial steps u^f = g < K and uq = g < K . □ 

As a result of Propositions 3.1, we can define 

{4 m = lim u l ' m , < I < L, < to < M. (3.16) 

n — >+oo 

It follows from Proposition 3.2 that u 1 ^ 1 < K. Letting n go to +oo, we can see from (3.10) that Uoo satisfies the 
difference equation (3.5). Therefore, 

Uoo = u. (3-17) 
In the following, we will study the convergence rate of {u n }n>o- 

Proposition 3.3. u n converges to u uniformly and 

max - u l r) < (1 - V A T (aTt)" K' (3 ' 18) 
where ij = rq^x+rlAt ^ * s a P os ^ ve constant. 

Proof. Let us dchnc 

e tsn = u i,m_~i mj E™ = m.axe l * n , E n = m.axE™. 

I m 

Proposition 3.1 and (3.17) ensure that e^ m is nonnegative. Moreover e n satisfies 



-p-e l -^ m + (1 +p ) e l r P + e l ^ m ~ ~ ^ ■ (^1) \ ^ 



(3.19) 



We can drop the first term on the left-hand-side of (3.19) because of the first inequality in (3.10) and e^ m being 
nonnegative. It gives us the inequality 

(1 +Po)e l A m [u l ' m - 9 l ] < [p~e l - 1 ' m + P+ e l + 1 ' m + e^ n+1 + XAtE n ^] [u l > m - g l ] , (3.20) 

in which we also used the assumption (3.9) to derive the upper bound for the integral term. 

If there are some (l,m) such that u ' n = g , since {u n }n>o is an increasing sequence from Proposition 3.1, we 
have u l,m = u l 1 i m for all n. Therefore, e 1 ^ 1 = for these (Z, m). On the other hand, if u> m > g l for some (/, m), we 
can divide u l,m — g l on both sides of (3.20) to get 

(1 +Po)e z /" < p-e l n - 1 ' m +p + e l + 1 ' m + e l n m+1 + XAtE^ 

(3.21) 

< p_E™ +p + E™ + + \AtE n _i. 

Since the right-hand-side of (3.21) does not depend on I, we can write 

K < VK +1 + (1 - ^y^-^-i, (3.22) 
A + r 

in which r\ = i + ( X +r)At e ■"•)• Note that (3.22) is also satisfied for all to, because even if u l ' m = g l for some 
(l,m), e\i m — as we proved above. If follows from (3.22) that 

K < v M-rn E M + (] _ _ ^ +„ + ... + j A (g ^ 

A + r 
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Since the terminal condition of u n , we have Eff = 0. Now maximizing the right-hand-side of (3.23) over to, we 
obtain that 

E n <(l- V M )-^E n - V 
X + r 

As a result, 

E„ < (1 -V M ) n (j-^-) ^o^O, asn^+oo. (3.24) 



□ 



Remark 3.2. As M 



i \ M 

1 - r) M = 1 - f \ , ) -+ 1 - e-^ T 7 



l + (A + r)T/M. 

which agree with the convergent rate (2.10) in the continuous case. 
Proposition 3.4. 

\u(x k ,mAt) - u{x k .mAt)\ -> 0, (3.25) 

as Ax, At, Az -> 0. 

Proof. Using the triangle inequality, let us write 

\u(x k ,mAt) - u(x k ,mAt)\ < \u(x k ,mAt) - u n (x k ,mAt)\ + \u n (x k ,mAt) - u n (x k ,mAt)\ 

+ \u n (x k , mAt) - u(x k , mAt) \ 

<K (l - e -(r+X)(T-rnAt } y \ " + „ . Q ^ + ^2 + (Az) 2) (3.26) 



K(l-r, M ) 



A 



A + 



for some positive constants K and K. The first and third terms on the right-hand-side of the second inequality are 
due to (2.10) and (3.18). The second term arises since the order of error from discretizing a PDE using implicit 
Euler scheme is 0((At) + (Ax) 2 ), the interpolation and discretization error from numerical integral are of order 
(Ax) 2 and (Az) 2 and the total error made at each step propagates at most linearly in n when we sequentially 
discretize (3.4). 

Letting At, Ax, Az — > in (3.26), we obtain that 

lim \u(x k ,mAt) - u(x k ,mAt)\ < ( K + k) ( -^-) ( 1 - e - (r+x)T ) " , 
At,Ax,Az^o \ /yA + r/V / 

in which we used (3.25). Since n is arbitrary the result follows. □ 

Remark 3.3. In Propositions 3.1 - 3.4, we have shown the convergence of the algorithm for completely implicit 
Euler scheme (6 = 1). In order to have the time discretization error as 0((At) 2 ), we will choose Crank- Nicols on 
scheme with 9=1/2 in the numerical experiments in the next section. From numerical results in Table 4, we shall 
see that Crank- Nicolson Scheme is also stable and the convergence is fast. 
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4. The Numerical Performance of the Proposed Numerical Algorithm 



In this section, we present the numerical performance of the algorithm proposed in the previous section. First, 
we compare the prices we obtain to the prices obtained in the literature. To demonstrate our competitiveness we 
also list the time it takes to obtain the prices for certain accuracy. We will use either the PSOR or the Brennan- 
Schwartz algorithm to solve the sparse linear system in (3.10); see Remark 4.1. All our computations are performed 
with C++ on a Pentium IV, 3.0 GHz machine. 

In Table 1, we take the jump distribution F to be the double exponential distribution 



We compare our performance with that of [17] and [16]. [17] obtain an approximate American option price formula, 
for by reducing the integro-pde equation V satisfied to a integro-ode following [5] . This approximation is accurate 
for small and large maturities. Also, they do not provide error bounds, the magnitude of which might depend 
on the parameters of the problem, therefore one might not be able to use this price approximation without the 
guidance of another numerical scheme. A more accurate numerical scheme using an approximation to the exercise 
boundary and Laplace transform was later developed by [16]. Our performance has the same order of magnitude 
as theirs. Our method's advantage is that it works for a more general jump distribution and we do not have to 
assume a double exponential distribution for jumps as [17] and [16] do. 

In Table 2 we compute the prices of American and European options in a Merton jump diffusion model, in which 
the jump distribution F is specified to be the Gaussian distribution 



We list the accuracy and time characteristics of the proposed numerical algorithm algorithm. We compare our 
prices to the ones obtained by [11, 12]. [11] used a penalty method to approximate the American option price, 
while we analyze the variational inequalities directly (see (3.5) and (3.10)). Moreover, our approximating sequence 
is monotone (see Proposition 3.1). 

In Table 3, We also list the approximated prices of Barrier options. We compare the prices we obtain with [19] 
where a Monte Carlo method is used. We do not list the time it takes for the alternative algorithms in Tables 2 
and 3 either because they are not listed in the original papers or they take unreasonably long time. 

In Table 4, we list the numerical convergence of the proposed algorithm algorithm with respect to grid sizes. We 
choose Crank-Nicolson scheme with 9 = 1/2 in (3.10) and solve the sparse linear system by either the Bernnan- 
Schwartz algorithm or the PSOR. 

Remark 4.1. Here we will analyze the complexity of our algorithm. Let us fix Ax/ At as a constant and choose the 
number of grid point in x to be N . For each time step, using the FFT to calculate the integral term in (3.10) costs 
0(N log N) computations. On the other hand, the Brennan- Schwartz algorithm, which uses the LU decomposition 
to solve the sparse linear system in (3.10) (see [15] pp. 283), needs 2N computations for each time step. However, 
PSOR needs C ■ N computations for each time step to solve (3.10) at each time step. Here, C is the number of 
iterations PSOR requires to converge to a fixed small error tolerance e. We will see in the following that PSOR is 
numerically more expensive than the Brennan- Schwartz algorithm. 




(4.1) 




(4.2) 
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For PSOR, the number of iterations C increases with respect to N . To see this, we start from the tri-diagonal 
matrix on the left-hand- side of (3.10) 



A 



1 + 6pa -6p+ 
-6p- 1 + 9 Po 



-8p+ 
-9p_ l + 8 Po J 

For the SOR (without projection) , the optimal relaxation parameter to is given by (see [22]) 

2 



where pj is the spectral radius of the Jacobi iteration matrix J = D^ 1 (A — D) with D as the diagonal matrix of A. 
Since p,j < || J||oo = @(p+ + + @Po), w e have 

lo < iv = 2 (4.3) 

i + Ji-IUIIi™ 



We will use ujq as the optimal relaxation parameter in our numerical experiments. On the other hand, since the 
largest eigenvalue X m ax of the SOR iteration matrix is bounded above by w— 1, using (3.6) and (4-3) we obtain that 



C = min{c > 0|(A roax ) c < e} = 0(VN). (4.4) 

Since 0(N 3 / 2 ) dominates O(NlogN), the complexity of the Iterated PSOR algorithm at each time step will be 
0(N^/ 2 ). Therefore, with O(N) time steps, the complexity for Iterated PSOR algorithm is 0(N b / 2 ). On the other 
hand, for the Iterated Brennan- Schwartz algorithm, since 0(N log N) dominates O(N), the complexity at each time 
step will be 0(N log N). Therefore, the complexity of the Iterated Brennan- Schwartz algorithm is 0(N 2 logiV) since 
we have 0(N) time steps. 

Please refer to Tables 1, 2, 3 and 4 for numerical performance of both algorithms. 

Next, we illustrate the behavior of the sequence of functions {v n {S,t)} n >Q and its limit V in Figures 1, 2 and 3. 
All the figures are obtained for an American put option in the case of the double exponential jump with K = 1 00, 
So = 100, T = 0.25, r = 0.05, a = 0.2, A = 3, p = 0.6, 771 = 25 and 772 = 25 (the same parameters are used in the 
8th row of Table 1) at a single run. 

Remark 4.2. (i) In Figure 1, we show, how V{S, 0) depends on the time to maturity, and that it fits smoothly 

to the put-pay-off function at s(0) (the exercise boundary). The y-axis is the difference between the option 
price and the pay-off function. As the time to maturity increases, the option price 1/(5,0) increases while 
the exercise boundary s(0) decreases. Even though the stock price process has jumps, the option price 
smoothly fits the pay-off function at s(0), as in the classical Black-Scholes case without the jumps. 
(ii) In Figure 2, we illustrate the convergence of the exercise boundaries t — > s n (t), n > 1. We can see from the 
figure that all s n (t) are convex functions. Also, the sequence {s„} ra >i is a monotone decreasing sequence, 
which implies that the continuation region is getting larger, and that the convergence of the free boundary 
sequence is fast. 

Moreover, we notice that, when the parameters are chosen such that (2.12) is satisfied, the free bound- 
aries are discontinuous at the maturity time. In addition, we have s(T—) = s n (T—) = S* < K, where S* 
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is the unique solution of (2.14). Furthermore, if F is the double exponential distribution as in (4-1), the 
integral equation (2.14) can be solved analytically to obtain 



With the parameters we choose, we get from (4-5) that S* = 98.39. It is close to our numerical result as 
one can see from Figure 2. 

(iii) In Figure 3, we illustrate the convergence of the sequence of prices {v n (S, 0)}„>o- Observe that this is a 
monotonically increasing sequence and it converges to its limit V{S, 0) very fast. 
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Table 1. Comparison between the proposed iterated jump algorithm with the method in [17] and 
[16], where the parameters were chosen as r = 0.05, S(0) = 100 and p = 0.6. Amin's price is 
calculated in [17] using the enhanced binomial tree method as in [3]. The accuracy of Amin's price 
is up to about a penny. The KPW 5EXP price from [16] is calculated on a Pentium IV, 1.8 GHz, 
while the iterated price is calculated on Pentium IV, 3.0GHz, both using C++ implementation. 
Run times are in seconds. For numerical algorithm we propose, the number of grid points in x is 
chosen as 2 6 and At = Ax. The option prices from both Iterated Brennan-Schwartz and Iterated 
PSOR are the same. Below "B-S" stands for the Brennan-Schwartz. 



American Put Double Exponential Jump Diffusion Model 





Parameter Values 




Amin's 


KW 


KPW 5EXP 




Proposed Algorithm 




K 


T 




\ 




72 


Price 


Value 


Error 


Value 


Error 


Time 


Value 


Error 


B-S Time PSOR Time 


90 


0.25 


0.2 


3 


25 


25 


0.75 


0.76 


0.01 


0.74 


-0.01 


3.21 


0.75 


o 


n ns 

u.uo 


0.12 


90 


0.25 


0.2 


3 


25 


50 


0.65 


0.66 


o 


0.65 


o 


3.25 


0.66 


0.01 


0.08 


0.12 


90 


0.25 


0.2 


3 


50 


25 


0.68 


0.69 


0.01 


0.68 


o 


2.97 


0.69 


0.01 


0.08 


0.12 


90 


0.25 


0.2 


3 


50 


50 


0.59 


0.60 


0.01 


0.59 





2.89 


0.59 





0.12 


0.12 


90 


0.25 


0.3 


3 


25 


25 


1.92 


1.93 


0.01 


1.92 





2.40 


1.93 


0.01 


0.09 


0.13 


90 


0.25 


0.2 


7 


25 


25 


1.03 


1.04 


0.01 


1.02 


-0.01 


3.18 


1.03 





0.12 


0.17 


90 


0.25 


0.3 


7 


25 


25 


2.19 


2.20 


0.01 


2.18 


-0.01 


2.97 


2.20 


0.01 


0.12 


0.20 


100 


0.25 


0.2 


3 


25 


25 


3.78 


3.78 





3.77 


-0.01 


3.08 


3.78 





0.12 


0.12 


100 


0.25 


0.2 


3 


25 


50 


3.66 


3.66 





3.65 


-0.01 


3.29 


3.66 





0.10 


0.12 


100 


0.25 


0.2 


3 


50 


25 


3.62 


3.62 





3.62 





2.88 


3.63 


0.01 


0.09 


0.12 


100 


0.25 


0.2 


3 


50 


50 


3.50 


3.50 





3.50 





3.00 


3.50 





0.13 


0.12 


100 


0.25 


0.3 


3 


25 


25 


5.63 


5.62 


-0.01 


5.63 





2.44 


5.63 





0.13 


0.15 


100 


0.25 


0.2 


7 


25 


25 


4.26 


4.27 


0.01 


4.26 





3.48 


4.27 


0.01 


0.17 


0.17 


100 


0.25 


0.3 


7 


25 


25 


5.99 


5.99 





5.99 





2.95 


6.00 


0.01 


0.17 


0.18 


90 


1 


0.2 


3 


25 


25 


2.91 


2.96 


0.05 


2.90 


-0.01 


2.43 


2.92 


-0.01 


0.63 


0.78 


90 


1 


0.2 


3 


25 


50 


2.70 


2.75 


0.05 


2.69 


-0.01 


2.38 


2.70 





0.69 


0.81 


90 


1 


0.2 


3 


50 


25 


2.66 


2.72 


0.06 


2.67 


0.01 


2.55 


2.68 


0.02 


0.64 


0.82 


90 


1 


0.2 


3 


50 


50 


2.46 


2.51 


0.05 


2.45 


-0.01 


2.30 


2.45 


-0.01 


0.68 


0.82 


90 


1 


0.3 


3 


25 


25 


5.79 


5.85 


0.06 


5.79 





2.48 


5.77 


-0.02 


0.70 


0.94 
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Table 2. Option price in Merton jump-diffusion model 

K=100, T=0.25, r=0.05, a — 0.15, A = 0.1. Stock price has lognormal jump distribution with fl = —0.9 and a = 0.45. 
For the iterated jump schemes, the number of grid points in x is chosen as 2 7 and At = Ax. Below "B-S" stands for the 
Brennan- Schwartz . 



Option Type a 


S(0) 


dFLV 6 




Proposed Algorithm 










Value 


Error 


LU(B-S) Time 


PSOR Time 


American Put 


90 


10.004 


10.004 c 





0.18 


0.24 




100 


3.241 


3.242 


0.001 








110 


1.420 


1.420 









European Put 


100 


3.149 


3.150 


0.001 


0.21 


0.18 


European Call 


90 


0.528 


0.528 





0.18 


0.18 




100 


4.391 


4.392 


0.001 








110 


12.643 


12.643 










"The option prices (for the same kind of option) for different S(0) are obtained from a single run. 
'The dFLV price comes from [11, 12]. 

c the option price is 10.001 using the iterated Brennan-Schwartz scheme. 

Table 3. European down-and-out barrier call option with Merton jump-diffusion model 

K=110, S(0)=100, T=l, r=0.05, a = 0.25, A = 2, rebate R=l, the Stock price has lognormal jump distribution with p, — 
and a = 0.1. For the algorithm we propose the number of grid points in x is chosen as 2 6 and At. = Ax. Below we use the 
acronyms LU or SOR to tell wheher we use the LU factorization or the SOR to solve for the sparse linear systems at each 
time step. 



Barrier H 


MA Price a 




Proposed Algorithm 






Value 


Error LU Time SOR Time 


85 


9.013 


8.988 


-0.025 0.52 0.71 


95 


5.303 


5.290 


-0.013 0.64 0.86 



"The MA price comes from [19] 
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Table 4. Convergence of the numerical algorithm with respect to grid sizes 

K=100, T=0.25, r=0.05, a = 0.15, A = 0.1, stock price has lognormal jump distribution with /2 = —0.9 and a = 0.45 (the 
same parameters that are used in [11]). The differential equation is discretized by the Crank-Nicolson scheme as (3.10) 
with 9 = 1/2. The logarithmic variable x — log S is equally spaced discretized on an interval [x m in,Xmax\ with Aa; = At. 
The numerical integral is truncated on the smallest interval [z m in, z max \, such that [x + jl — 4<r, x + u + 4er] will be inside 
[zmin, Zmax] for any x £ [xmin, Xmax]- The step length for the numerical integral is chosen the same as the step length in x, 
i.e. Az = Aa-. The number of grid points for to implement the FFT is chosen as an integral power of 2. The error tolerance 
for PSOR method is 10 -8 and for the global iteration is 10 -6 . Run times are in seconds. Each row in the "Difference" 
column of the following table is vpsor(L, M) — vpsor(L/2, M/2). "B-S" stands for the Brennan-Schwartz algorithm. The 
number of global iteration is 3 for all the following numerical experiments. 



S(0) 


No. of grid 


No. of time 


B-S Value 


B-S 


PSOR Value 


Difference 




PSOR 


Max. No. of 




points in x ( L ) 


steps ( M ) 


VB-S 


Time 


VPSOR 






Time 


PSOR iterations 




64 


30 


10.00230 


0.06 


10.00573 


n.a. 




0.06 


16 




128 


58 


10.00142 


0.21 


10.00429 


-0.00144 




0.24 


21 


90 


256 


115 


10.00192 


0.84 


10.00396 


-0.00033 




0.99 


28 




512 


230 


10.00218 


3.51 


10.00387 


-0.00009 




4.50 


39 




64 


30 


3.24074 


0.06 


3.24465 


n.a. 




0.06 


16 




128 


58 


3.24008 


0.21 


3.24180 


-0.00285 




0.24 


21 


100 






















256 


115 


3.24046 


0.84 


3.24115 


-0.00065 




0.99 


28 




512 


230 


3.24058 


3.51 


3.24103 


-0.00012 




4.50 


39 




64 


30 


1.42048 


0.06 


1.42146 


n.a. 




0.06 


16 




128 


58 


1.41941 


0.21 


1.41991 


-0.00155 




0.24 


21 


110 






















256 


115 


1.41958 


0.84 


1.41966 


-0.00025 




0.99 


28 




512 


230 


1.41962 


3.51 


1.41960 


-0.00006 




4.50 


39 



Using (4.4), the number of SOR iterations can be calculated. The calculation gives 11, 16, 22 and 31. Comparing with the 
last column of above table, the maximum numbers of PSOR iteration are slightly larger than these theoretical predicted 
SOR iteration times. Moreover, when L = 512 the the ratio between the maximum number of PSOR iteration and VX is 
1.72. This confirms the analysis in Remark 4.1 that the maximal PSOR iteration time grows as the order of s[L. 
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The parameters for the following three figures are K = 100, So = 100, T = 0.25, r = 0.05, a = 0.2, A = 3, the 
stock price has double exponential jump with p = 0.6, r/i = 25 and r/2 = 25 (the same parameters used in the 8th 

row of Table 1). 

Figure 1. The option price function S — ► V(S, 0) smoothly fits the pay-off function (K — S) + 
at s(0). V(S, 0) increases and s(0) (V(S, 0) — (K — S) + = at s(0)) decreases as time to maturity 
T increases. 



T: time to 
maturity 

T = 0.25 




Stock Price ( S ) 
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Figure 2. Iteration of the Exercise Boundary: s„ (t) J, s(t), t e [0, T). Both s n {t) and s(i) will 
converge to S* < K as t — > T. 




0.05 0.1 0.15 0.2 0.25 
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